1. Field of the Invention
The present invention relates to a simulation system which performs a numerical simulation of a system that suffers random disturbance.
2. Description of the Related Art
In the field of the natural sciences, it is generally difficult to analytically solve a mathematical model that describes a physical system. In most cases, the initial-value problem in the mathematical model is solved approximately by means of numerical processing. With the recent advance of information processing technology, a wide variety of simulators adapted to perform such numerical processing have been devised. However, these simulators require an overwhelming amount of computation so as to increase the accuracy of the simulation, which may decrease simulation versatility. Therefore, the important question is how to build reliable and versatile simulators.
In the field of engineering which deals with the behavior of a physical system such as an industrial plant or a robot, it is an important problem to estimate the behavior of the physical system and control it on the basis of observed data from sensors. This estimation is generally performed by a numerical simulation using a simulator.
The performance of a simulator affects the reliability of the result of a numerical simulation, which, in turn, determines the performance of a controller which controls the physical system. For this reason, various simulators have been developed so far according to the type of a physical system which is an object of numerical simulation and a reliability required of a controller.
A model of a physical system to be processed is derived on the basis of the laws of nature. In this case, the effect of random disturbance such as noise in many actual physical systems is often neglected so that they are idealized or simplified. Depending on the physical system, the effect of random disturbance sometimes cannot be neglected.
In designing a controller for an artificial satellite by way of example, it is impossible to neglect random disturbance such as radiation pressure of the sun, gravity, or collision with the wreck of another satellite drifting in space. With a controller for an industrial plant or robot, it is sometimes impossible to neglect random disturbance such as noise mixed in signals from sensors for observing a controlled object.
In a simulation that aims at analysis and control of a physical system which undergoes the effect of random disturbance, it is required to provide a physical model that reflects the random disturbance and to process it numerically.
In general, the model of any physical system can be represented by the ordinary differential equation: EQU x=f(x,t), x(t.sub.0)=x.sub.0 ( 1)
where x=x(t) indicates a physical quantity for simulation, such as displacement, temperature or the like, in the physical system and varies with time t. x indicates differentiation of x with respect to t. f(x, t) indicates a function of x and t, and x.sub.0 is the initial value of x at time t=t.sub.0.
In a numerical simulation for the physical model represented by equation (1), let the approximate value of physical quantity x(t.sub.n) at t.sub.n =t.sub.0 +nh be x.sub.n. Then, the approximate value x.sub.n+1 of the physical quantity x(t.sub.n+1) at t.sub.n+1 will be calculated by the equation: EQU x.sub.n+1 =x.sub.n +.PSI.(x.sub.n, t.sub.n, h) (2)
Note that x.sub.0 is not an approximate value but the initial value of x at time t=t.sub.0. Here, h indicates a timing step size or an internal of the simulation and .PSI.(x.sub.n, t.sub.n, h) provides an increment (or a time increment) of the approximate value x.sub.n at t.sub.n. How .PSI.(x.sub.n, t.sub.n, h) is to be set up with x.sub.n, t.sub.n and h is so important as to affect the reliability of the result of numerical simulation.
Methods of setting up .PSI.(x.sub.n, t.sub.n, h) in equation (2) can be roughly classified into two methods: Taylor series method and Runge-Kutta method. The Taylor series method, which is in principle simple, provides .PSI.(x.sub.n, t.sub.n, h) by Taylor series terms up to the order required for simulation. With this method, however, the coefficient of the p-th order term generally will involve the p-th order derivative .differential..sup.p f(x, t)/.differential.x.sup.p of the function f(x, t). If the function f(x, t) is of a simple type, for example, a linear function of x, then that derivative will be obtained easily. If, however, f(x, t) is a complex non-linear function, then it will be difficult to provide its high-order derivative as a function of x and t. With the Taylor series method, it is required for an operator to find derivatives of f(x, t) and enter them into a simulator. Furthermore, the computation of derivatives is a heavy burden on a simulator. Thus, the Taylor series method is poor in versatility. In numerical simulation, therefore, the Runge-Kutta method is often used, which forms .PSI.(x.sub.n, t.sub.n, h) by the equation: ##EQU1## where s is a natural number, and k.sub.i (k parameter) is a value of function f(x, t) which corresponds to appropriately chosen x and t. The k parameter k.sub.i is generated by the following computational processing. ##EQU2##
According to equations (4), k.sub.i can be obtained explicitly using k.sub.1 to k.sub.i-1. Thus, this method is called an explicit Runge-Kutta method. In equations (3) and (4), the b parameter b.sub.i and the a parameter a.sub.ij are parameters indicating weights for the k parameters, which are stored in memories within a simulator. Depending on how b parameters and a parameters are chosen, the incremental value .PSI.(x.sub.n, t.sub.n, h) for the approximate value x.sub.n varies, which affects the accuracy of the approximate value x.sub.n+1 and consequently determines the performance of the simulator.
With the method using k parameters, there is no need of providing any differential form of the coefficient function f(x, t) in performing the numerical processing of equation (2). Therefore, this method is far more versatile than the Taylor series method.
Assuming the solution of equation (1) at time t=tn to be x(t.sub.n), a difference or error between the approximate value x.sub.n based on numerical simulation and x(t.sub.n) is given by x(t.sub.n)-x.sub.n. The smaller the error is, the higher the reliability of the results of simulation becomes. The timing step size h serves as a measure for estimating the magnitude of that error. When it is assured that the error is on the order of h.sup.r+1, the simulator or the result of simulation is said to have an approximation accuracy of r order. The greater the value of r is, the higher the reliability of simulation becomes.
When a physical system to be simulated undergoes the effect of random disturbance such as noise, the physical model for that system is described using the following stochastic differential equation of Ito's type instead of equation (1). EQU x(t)=f(x)+g(x)w(t), x(t.sub.0)=x.sub.0 ( 5)
where x=x(t) stands for the state quantity in the physical system to be simulated, f(x) and g(x) are each a function of x, and w(t) represents a noise model the value of which depends on the stochastic process of a Brownian motion model. This model includes a second term associated with random disturbance added to the right-hand side of equation (1). x.sub.0 is the initial value of x at time t=t.sub.0.
Since equation (5), representing a physical model subjected to random disturbance, unlike equation (1), has an added term associated with random disturbance, .PSI.(x.sub.n, t.sub.n, h) in equations (3) and (4) cannot be used as it is. In order to handle the stochastic-process-dependent random disturbance term, it is required to expand .PSI.(x.sub.n, t.sub.n, h) in equations (3) and (4) to the case of a stochastic differential equation. Conventional methods of setting up the expanded an increment .PSI.(x.sub.n, t.sub.n, h) using k parameters include a stochastic Heun method, a 3-stage explicit method that assures 1.5 order accuracy, and an asymptotically efficient method.
The stochastic Heun method (RUMELIN, SIAM J. Vol.19, No.3, June 1982) uses a stochastic-process-dependent term added to the right-hand side of equation (3) as .PSI.(x.sub.n, t.sub.n, h) for the physical model described by equation (5). This method is a direct expansion of equation (3) for a physical model subjected to stochastic-process-dependent random disturbance, assuring first order approximation accuracy. Unlike equation (3), this expansion is characterized in that the k parameters include a value of function g(x) and a value of a first derivative of function g(x).
The 3-stage explicit method (SAITO, MITSUI, Japan applied mathematics society journal Vol. 2, No. 1, 1992) further expands .PSI.(x.sub.n, t.sub.n, h) in the stochastic Heun method using two types of stochastic processes. In general, this method assures a high approximation accuracy of 1.5 order, and the k parameters include a value of a first derivative of each of functions f(x) and g(x) and a value of a second derivative of g(x).
The asymptotically efficient method (NEWTON, SIAM J. Vol. 51, No. 2, April 1991) further adds a term relating to (h).sup.1/2 to .PSI.(x.sub.n, t.sub.n, h), thereby simplifying computational processing of the k parameters. With this method, the k parameters include no value of first derivative of function g(x) and can be calculated by using values of f(x) and g(x). For this reason, this method is superior in simulation versatility. The approximation accuracy of this method is one order as in the case with the stochastic Heun method.
The numerical simulation based on the conventional time-increment forming methods have the following problems.
With the stochastic Heun method, in calculating the increment .PSI.(x.sub.n, t.sub.n, h) for the approximate value x.sub.n, two sets of k parameters are calculated using the x.sub.n value already obtained, t.sub.n and h, and the k parameters are weighted suitably and then added. At this point, a calculation expression for a set of k parameters includes a first derivative of the function g(x) describing random disturbance. In performing a numerical simulation, therefore, the value of the first derivative of the function g(x) must be calculated. For this reason, as with the Taylor series method, the amount of calculation increases, increasing the load on the simulator. Thus, this method is not versatile. In addition, the approximation accuracy of the simulation remains at first order and hence the reliability of the simulation is not necessarily high.
The 3-stage explicit method is more reliable than the stochastic Heun method for simulation results. However, since the k-parameter calculating expression includes first derivatives of the functions f(x) and g(x) and a second derivative of g(x), the function forms of those derivatives must be given. The necessity of not only the values of the first derivatives of f(x) and g(x) but also the value of the second derivative of g(x) results in a significant increase in the amount of calculation in comparison with the stochastic Heun method. In practice, it is difficult to apply this method to simulators.
In comparison with the stochastic Heun method and the 3-stage explicit method, the asymptotically efficient method has an advantage that it does not require provision of derivatives of g(x) and the like in performing a numerical simulation. With this method, the k-parameter calculation expression is described by the functions f(x) and g(x). Thus, this method permits numerical processing to be performed easily and is superior in versatility. However, the drawback of this method is low reliability because the approximation accuracy of simulation remains at first order.
To simulate a system subjected to the random-disturbance effect, it is essential to develop a simulator that finds numerically the solution of the stochastic differential equation described by equation (5). Nevertheless, few practical simulators have been developed so far because the conventional simulation methods have the above-described problems with respect to reliability and versatility.